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1. INTRODUCTION 

Precipitation kinetics in alloys, like spinodal decomposition, nucleation and growth, 
or phase ordering, are now often studied at an atomistic scale using Monte Carlo 
simulations. So as to be able to reproduce the different kinetic behaviors during these 
transformations, one needs to adopt a realistic description of the diffusion. Therefore it 
is better to use a vacancy-diffusion mechanism than a direct atom exchange mechanism. 
It is then possible to explain why different kinetic pathways are observed. For instance, 
the vacancy diffusion mechanism can predict the importance of monomer diffusion with 
respect to the diffusion of small clusters [1-4]. This leads to a difference in the cluster 
size distribution during precipitation [2] and determines the coarsening mechanism 
(evaporation-condensation or coagulation) [3,4]. One can predict too the slowdown of 
precipitation kinetics by vacancy trapping due to the addition of a third component 
impurity [1]. Finally, different interactions of solute atoms with vacancy lead to a 
difference of precipitate / matrix interface morphology during the kinetic pathway, the 
interface being diffuse for a repulsion between vacancy and solute atoms and sharp for 
an attraction [4]. 

One drawback of kinetic Monte Carlo simulations using vacancy-diffusion mecha- 
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nism is that they hmit themselves to pair interactions to describe configurational energy 
of alloys. Multisite interactions including more than two lattice sites are necessary if 
one wants to fully reproduce the thermodynamics of a system [5,6]. These interactions 
reflect dependence of bonds with their local environment and as a consequence break 
the symmetry imposed by pair interactions on phase diagram. It is interesting to notice 
that in Calphad approach [7] one naturally considers such interactions by describing for- 
mation energy of solid solutions with Redlich-Kister polynomials, and that coefficients 
of these polynomials can be mapped onto an Ising model to give effective interactions 
including more than two lattice sites [8,9]. Moreover these interactions allow one to 
understand shapes of precipitates in alloys [10] and can lead to a prediction of coher- 
ent interface energy in really good agreement with ab-initio calculations performed on 
supercells [11]. Nevertheless, to study kinetics in Monte Carlo simulations with such 
interactions, one usually uses a direct atom exchange mechanism [12], and thus looses 
all kinetic effects due to vacancy-diffusion mechanism. 

We incorporate these multisite interactions in a kinetic model using a vacancy- 
diffusion mechanism to study precipitation kinetics of AlaZr in Al-Zr solid solution. For 
small supersaturation in zirconium of the aluminum solid solution, it has been experi- 
mentally observed that AlaZr precipitates are in the metastable LI2 structure [13-15]. 
These precipitates are found to have mainly spherical shape (diameter ~ 10 — 20 nm), 
as well as rod-like shape [13]. For supersaturation higher than the peritectic concentra- 
tion, nucleation is homogeneous and precipitates are coherent with the matrix [13,14]. 
With prolonged heat treatment, if the temperature is high enough, the metastable LI2 
structure can transform to the stable one DO23. Using a phase field method, Proville 
and Finel [16] modelled these two steps of the precipitation, i.e. the transient nucle- 
ation of the LI2 structure and the transformation to the DO23 structure. In this work, 
we only focus on the precipitation first stage, where AlaZr precipitates have the LI2 
structure and are coherent with the matrix. 

We first use ab-initio calculations to fit a generalized Ising model describing thermo- 
dynamics of Al-Zr system. We then extend description of the configurational energy of 
the binary Al-Zr to the one of the ternary Al-Zr- Vacancy system and adopt a vacancy- 
atom exchange mechanism to describe kinetics. This atomistic model is used in Monte 
Carlo simulations to study diffusion in the Al-Zr solid solution as well as precipitation 
kinetics of AlsZr. We mainly focus our study on detecting any influence of multisite 
interactions on kinetics. 

2. THERMODYNAMICS OF Al-Zr BINARY 

2.1. Ab inito calculations 

We use the full-potential linear- muffin-tin-orbital (FP-LMTO) method [17-19] to 
calculate formation energies of different compounds in the Al-Zr binary system, all 
based on a fee lattice. Details of ab initio calculation can be found in appendix^ They 
are the same as in our previous work [20], except the fact that we use the generalized 
gradient approximation (GGA) instead of the local density approximation (LDA) for 
the exchange- correlation functional. 

The use of GGA for the exchange correlation energy leads to a slightly better 
description of the Al-Zr system. The approximation does not fail to predict phase 
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stability of pure Zr [21] as LDA does: if one does not include generalized-gradient 
corrections, the stable structure at K for Zr is found to be the uj one (hexagonal with 
3 atoms per unit cell) and not the hep structure. 

Another change depending on the approximation used for the exchange-correlation 
functional is that formation enthalpies obtained with GGA for the different Al-Zr com- 
pounds are a little bit lower (a few percent) than with LDA. For the DO23 structure of 
AlsZr (table P), GGA predicts a formation energy which perfectly reproduces the one 
measured by calorimetry [22]: including generalized-gradient corrections has improved 
the agreement. The energy of transformation from the LI2 to the DO23 structure, 
AE = —23 meV/atom, agrees really well too with the experimental one measured 
by Desh et al. [23], but this was already true with LDA. Gradient corrections have 
improved the agreement for the equilibrium volumes too: with the LDA, they were 
too low compared to the available experimental ones. Considering the values of the 
relaxed equilibrium parameters, shape of the unit cell and atomic positions, no change 
is observed according to the approximation used, both LDA and GGA being in good 
agreement with measured parameters. 

Table 1: Calculated equilibrium volumes Vq, c'/a ratios (c' = c/2 for the DO22 phase 
and c' = c/4 for the DO23 phase), atomic displacements (normalized by a), and energies 
of formation for Al3Zr compared to experimental data. 







Vo 

(A^/atom) 


c'/a 


Atomic 
displacements 


AE 

(eV/atom) 


L12 


GGA" 


16.89 






-0.478 




LDA^ 


16.12 






-0.524 


D022 


GGA" 


17.40 


1.138 




-0.471 




LDA^ 


16.60 


1.141 




-0.525 


D023 


GGA" 


17.16 


1.080 


5ai = +0.0013 


-0.502 










= -0.0239 






LDA* 


16.35 


1.087 


5ai = -0.0021 


-0.548 










= -0.0273 






Exp.'^ 


17.25 


1.0775 


= +0.0004 












5zr = -0.0272 






Exp.-^ 








-0.502 + 0.014 



TP-LMTO calculations (present work) 
''FP-LMTO calculations [20] 
■^Neutron diffraction [24] 
''Calorimetry [22] 



2.2. Cluster expansion of the formation energy 

In order to express the formation energy of any Al-Zr compound based on a perfect 
fee lattice, we make a cluster expansion [25] of our FP-LMTO calculations to fit a 
generalized Ising model. This allows us to obtain the energy of any configuration of 
the fee lattice. 
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Considering a binary alloy of N sites on a rigid lattice, its configuration can be 
described through an Ising model by the vector cr = {ai, (72, ... , ajy} where the pseudo- 
spin configuration variable cxj is equal to ±1 if an A or B atom occupies the site i. Any 
structure is then defined by its density matrix p**, p'*(cr) being the probability of finding 
the structure s in the configuration cr. 

With any cluster of n lattice points a = {ii,i2, ■ ■ ■ yin} we associate the multisite 
correlation function 



where the sum has to be performed over the 2 possible configurations of the lattice. 

Clusters related by a translation or a symmetry operation of the point group of the 
structure have the same correlation functions. Denoting by the number of such 
equivalent clusters per lattice site, or degeneracy, the energy, like any other configura- 
tional function, can be expanded in the form [25] 



where the sum has to be performed over all non equivalent clusters and the cluster 
interaction Ja is independent of the structure. 

The cluster expansion as defined by equation |21 cannot be used directly: a truncated 
approximation of the sum has to be used. The truncation is made with respect to the 
number of points contained in a cluster, thus assuming that order effects on energy are 
limited to a small set of lattice points. It is truncated too with respect to distance 
between sites. Long range interactions are important mostly if one wants to fully 
reproduce elastic effects [26]. 

We use in the expansion of the energy six different clusters: the empty cluster {0}, 
the point cluster {!}, the pairs of first and second nearest neighbors {2,1} and {2,2}, the 
triangle of first nearest neighbors {3,1}, and the tetrahedron of first nearest neighbors 
{4,1}. The corresponding cluster interactions are obtained by making a least square fit 
of compound energies calculated with FP-LMTO. All 17 used compounds are lying on a 
perfect fee lattice: energies are calculated without any relaxation of the volume, of the 
shape of the unit cell, or of the atomic positions. The lattice parameter used is the one 
which minimizes the cohesive energy of pure Al, a = qai = 4.044 A. We choose to fit the 
cluster expansion for the equilibrium lattice parameter of Al because we are interested in 
describing thermodynamics of the Al rich solid solution as well as of AlsZr precipitates 
in the LI2 structure. These precipitates have an equilibrium lattice parameter close to 
the one of pure Al, a = 4.073 A as obtained from FP-LMTO calculations with GGA, 
and during the nucleation stage they are coherent with the Al matrix. Consequently, 
such an expansion should be able to give a reasonable thermodynamic description of 
the different configurations reached during this precipitation stage where precipitates 
are coherent. 

Coefficients of the cluster expansion of the energy are given in table El Comparing 
the values of the many-body interactions, we see that the main contribution to the 
energy arises from the pair interactions and that the 3- and 4-point cluster contributions 
are only corrections. Signs of pair interactions reflect the tendency of Al and Zr atoms to 




(1) 




(2) 



a 
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Table 2: Formation energy relative to pure fee elements for Al-Zr compounds lying on 
a perfect fee lattice (a = qai = 4.044 A) obtained from a direct FP-LMTO calculations 
and from its cluster expansion. 





Pearson 


Structure 


^/orm (eV/atom) 




symbol 


type 


FP-LMTO 


CE 


Al (fee) 


cF4 


Cu 


0. 


0. 


Al4Zr (DIJ 


tllO 


MoNi4 


-0.421 


-0.491 


AlgZr (LI2) 


cP4 


CuaAu 


-0.728 


-0.671 


AlgZr (DO22) 


tI8 


AlgTi 


-0.617 


-0.643 


AlaZr (DO23) 


tI16 


AlgZr 


-0.690 


-0.657 


Al2Zr iP) 


tI6 


MoSia 


-0.513 


-0.482 


AlZr (Llo) 


tP4 


AuCu 


-0.803 


-0.780 


AlZr (Lli) 


liR32 


CuPt 


-0.448 


-0.466 


AlZr (CH40) 


tis 


NbP 


-0.643 


-0.723 


AlZr (D4) 


cF32 


7a 


-0.489 


-0.414 


AlZr (Z2) 


tP8 


7a 


-0.345 


-0.333 


Zr2Al 


tI6 


MoSi2 


-0.443 


-0.482 


ZrgAl (LI2) 


cP4 


CusAu 


-0.640 


-0.603 


ZraAl (DO22) 


tI8 


AlgTi 


-0.570 


-0.574 


ZrsAl (DO23) 


tI16 


AlgZr 


-0.603 


-0.589 


Zr4Al (DIJ 


tllO 


MoNi4 


-0.390 


-0.437 


Zr (fee) 


cF4 


Cu 


0. 


0. 



Description of structures D4 and Z2 can be found in Ref. [6] 

Table 3: Cluster expansion of the formation energy. 



Cluster 


Da 


J a 

(eV/atom) 




{0} 


1 


-4.853 




{1} 


1 


0.933 




{2,1} 


6 


97.5 X 10" 


-3 


{2,2} 


3 


-28.4 X 10" 


-3 


{3,1} 


8 


4.2 X 10- 


-3 


{4,1} 


2 


13.1 X 10" 


-3 



form heteroatomic first nearest neighbor pairs and homoatomic second nearest neighbor 
pairs. 

In tableEl we compare the formation energies of the different compounds directly ob- 
tained from FP-LMTO calculations with the ones given by their cluster expansion. The 
standard deviation equals 41 meV/atom and the maximal difference is 79 meV/atom. 
This could have been improved by including more clusters in the expansion of the en- 
ergy or by using a mixed-space cluster expansion [26]. Nevertheless, this would not 
have changed the main characteristics of the Al-Zr system, i.e. the short range order 
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tendency given by pair interactions, as well as the dependence on local environment of 
the interactions given by 3- and 4-point cluster interactions. In order to be able to build 
a realistic kinetic model and to run Monte Carlo simulations in a reasonable amount 
of time, we have to keep the thermodynamic description of Al-Zr system as simple as 
it can be. Therefore we do not try to improve expansion convergence and we focus our 
work on the influence of the 3- and 4-point cluster interactions on the thermodynamic 
and kinetic properties. 

2.3. Phase diagram 

T (K) I \ \ \ \ \ \ \ \ \ 




10 20 30 

%Zr 



Fi gure 1: Al rich part of the phase diagram corresponding to the equihbrium between the fee sohd 
solution and the LI2 structure given by our set of parameters (table OJ- (a) Comparison of the phase 
diagrams obtained with pair, triangle, and tetrahedron interactions (solid line) and the one obtained 
with only pair interactions (dotted line), (b) Comparison with the predicted metastable solubility 
limit [20] (dashed line). 

We use the cluster-variation method (CVM) [27] in the tetrahedron-octahedron 
(TO) approximation [28, 29] to study the equilibrium between the fee Al-rich solid 
solution and the LI2 structure (Fig. ^ corresponding to energy parameters of table El 
At low temperature, the 2 sublattices of the LI2 structure remain highly ordered, as 
at the experimental peritectic melting temperature (T ~ 934 K) the Zr concentrations 
of the two sublattices are respectively 100. and 1.8 at.%. Turning out the energy 
coefficients of the first nearest neighbor triangle and tetrahedron (J3 = J4 = 0), we 
see that these many-body interactions have a thermodynamic infiuence only at high 
temperature (Fig.H(a)), as for temperatures below 1000 K the phase diagram remains 
unchanged with or without these interactions. 

In figure n (b), we compare the Zr solubility limit in the fee solid solution corre- 
sponding to the present work energy parameters with our previous estimation of this 
metastable solubility limit [20]. We should point out that the solubility limit obtained 
with the parameters given by table El corresponds to a coherent equilibrium between the 
fee solid solution and the LI2 structure as the energy coefficients of the expansion have 
been calculated for a perfect fee lattice at the parameter of pure Al. This leads to a 
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destabilization of the ordered phase and this is the main reason why we obtain a higher 
solubihty than the estimated one corresponding to the equihbrium between incoherent 
phases. Another reason is that we use the cluster expansion to compute AlaZr cohesive 
energy, and thus get a small error on this energy, whereas in our previous study we 
directly used the value given by FP-LMTO calculation. 

3. KINETIC MODEL 

In order to be able to build an atomistic kinetic model, we have to generalize our 
thermodynamic description of the Al-Zr binary system to the one of the Al-Zr- Vacancy 
ternary system. To do so, we recast first the spin-like formalism of the cluster expansion 
into the more convenient one of the lattice gas formulation using occupation numbers 
[5] . This will allow us to obtain effective interactions for the different configurations of 
the tetrahedron of first nearest neighbors and of the pair of second nearest neighbors. 
Atom-vacancy interactions can then be introduced quite easily. 

3.1. Effective interactions 

Instead of using the pseudo-spin variables (T„ as we did in chap 12. 2. [ this will be 
easier for the following to work with occupation numbers p^, being equal to 1 if an 
atom of type i occupies the site n and to otherwise. In a binary alloy, occupation 
numbers and pseudo-spin variable at site n are related by 

Pn = and = ^ (3) 

For the Al-Zr binary system, we included in our truncated cluster expansion of 
the energy first nearest neighbor interactions up to the pair, triangle, and tetrahedron 
clusters and a second nearest neighbor pair interaction. Thus, using the occupation 
numbers p^, the expression of the energy becomes 

^ n,m,p,q ^ r,s 

i,j,k,l i,j 

where the first sum runs over all sites {n,m,p,q) forming a first nearest neighbor 
tetrahedron and all their different configurations k, I), and the second sum over all 
sites (r, s) forming a second nearest neighbor pair and all their different configurations 
Ng is the number of lattice sites, e-j^^ the effective energy of a first nearest 

(2) 

neighbor tetrahedron in the configuration {i,j,k,l), and e^y the effective energy of a 
second nearest neighbor pair in the configuration 

Writing the energy with these effective interactions increases the number of depen- 
dent variables. Therefore several choices of these effective energies correspond to the 
same cluster expansion, then to the same thermodynamic and kinetic properties. If we 
make the assumption that second nearest neighbor interactions do not contribute to the 

(2) (2) 

cohesive energy of pure elements, i.e. e\-^ = and e^^ = 0, we obtain as many effec- 
tive interactions as parameters in the truncated cluster expansion. Such an assumption 
does not have any physical influence and it just guarantees that homo-atomic effective 
interactions, e^]^^^ and e^^\, do not depend on the on the nature of B atom. Effective 
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energies of the first nearest neighbor tetrahedron in its different configurations are then 
related to the cluster expansion coefficients by the equations 



/ ^AAAA \ 
^AAAB 
^AABB 
''ABBB 
\ ^BBBB 



1 

12 



/6 


6 


6 


6 


6\ 


/ -Do Jo + -02,2^2,2 


\ 


6 


3 





-3 


-6 








6 





-2 





6 




D2,lJ2,l 




6 


-3 





3 


-6 




-D3,1<^3,1 




V6 


-6 


6 


-6 


6/ 


v 


-04,1 J4,i 


/ 



(5) 



and the second nearest neighbor pair interaction by the equation 



.(2) 
-AB 



"2-02,2^2,2- 



(6) 



For Al-Zr binary system, tetrahedron effective interactions corresponding to the cluster 
expansion of chap. ITD can be found in tabled two sets are given depending if J3 and 
J4 are taken from the cluster expansion of table 01 or are supposed equal to zero. For 



both sets e 



(2) 
AB 



+0.057 eV. 



3.2. Decomposition of effective interactions 

As we wrote before, several sets of effective interactions produce the same cluster 
expansion. In the following, we generate the set of interactions useful for our kinetic 
model for which we have to count bonds we break for vacancy-atom exchange. 

Different contributions are included in the effective energy e^^j],;. One part of the 
energy is due to the bonding corresponding to the six different pairs of atoms contained 
in the tetrahedron, each of these pairs belonging to two different tetrahedrons. Then one 
has to add corrections due to order on the four triangles contained in the tetrahedron 
and another correction due to order on the tetrahedron itself. This decomposition leads 
to the relation 

^(1) _ 1 , ,(1) , ,(1) , ,(1) , ,(1) , 

where e^-j^ is the effective energy of the ffist nearest neighbor pair in the configuration 

(i, j) and ejj|, and e^j^^ the corrections to add to pair energy due to order on triangles 
and on the tetrahedron 

Using the previous breakdown of the tetrahedron effective energy, the expression |3] 
of the energy becomes 

~ TTTT ^ij^PnPm + TTTT '^ijkPnPmPp 

■^^^s n,m "^^^s n,m,p 

i,j i,j,k 

+ W ^ ^llPnP'-P'A + ^ E ^P^Pi- 
* n,m,p,q * r,s 
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As this is just another mathematical way to rewrite the cluster expansion |21 of the 
energy, the following relations holds : 

3 ^,(1) , o.(i) , .(1) 



D,J, + D2,2J2,2 = 2 + 2e^B + 4b j (9a) 

^^AAA + "^^AAB + '^^ABB + ^BBB 

Ah^) +4?(i) +6?^'^ +4?^') +^^^ 
+ g y^AAAA + ^^AAAB + ^^AABB + ^^ABBB + ^BBBB 

DiJi = Sfei^i-egi) (9b) 






AAB ^ABB ^BBB 
"T 2 \^ AAAA "1" ■^'^ABBB '^BB 

4 V^^^^ ^*=AA_BB ^ '^BBBB 
-03,1^3,1 = e^^AA - 3e^Ai? + ^^ABB - ^BBB (9d) 

+ _2?^^) +2?^') ' 

2 V ^^"^^ ^^^^ ^^^^ BBBB ^ 

n T _ -47(1) ^ (-Qp^ 

J^4,l-'4,l — g ^'^AAAB ~^ ^'^AABB ^^ABBB^^BBBBJ W^J 

As we want e|j|, to be the energetic corrections due to order on triangles, the con- 
tribution to J2,i of e^]j^, ^AABi ■ ■ ■ must equal zero (second term in right hand side 
of eq. IHEl). For the same reason, the contribution to J2,i and J^^i of '^aaaai '^aaabi 
must equal zero (last term in right hand side of eq. Ec] and l9dj) . We require 
too that triangle and tetrahedron order corrections do not contribute to the cohesive 
energy of pure elements, as we did for second nearest neighbor pair interactions. Thus, 
'^AAA ~ '^BBB ~ and '^AAAA ~ '^BBBB ~ ^- With these restrictions, all parameters 
entering in the expression |H1 of the energy are well determined. 

The first nearest neighbor pair effective energies are thus 

^AA = ^(^0^0 + ^2,1^2,1 + ^2,2^2,2 + ^3,1^3,1 + 1^4,1^4,1) (10a) 

i {DqJq - ^2,1 J2,l + ^2,2^2,2 + A,1^4,l) (10b) 
D 

i (DoJo - DiJi + ^2,1 J2,l + ^^2,2^2,2 " ^^3,1 ^3,1 + ^^4,1 ^4,0 , (lOc) 
D 

the order corrections on first nearest neighbor triangle 

- - / nn 

^AAB - ~^ABB - -g-'^S.l-'S,!, IJ-J-j 

and the order corrections on first nearest neighbor tetrahedron 

'^AAAB = ^ABBB = — -D4_i J4_i (12a) 
^AABB = 0. (12b) 



^AB 

,(1) 
-BB 
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Inverting the system El one can easily express all these quantities from the effective 
tetrahedron energies eijki too. 

3.3. Interactions with vacancy 

Within the previous formalism, we can easily introduce atom-vacancy interactions. 
These interactions are a simple way to take into account the electronic relaxations 
around the vacancy. Without them, the vacancy formation energy Ey"^ in a pure metal 
would necessarily equal the cohesive energy = 0.69 eV [30] and E""'' = 3.36 eV 

for fee Al). 

We only consider first-nearest neighbor interactions with vacancies and we do not 
include any order correction on triangle and tetrahedron configurations containing at 
least one vacancy, i.e. e^jy = e|j],y = where i, j, and k are any of the species Al, Zr, 
and V. The vacancy formation energy in a pure metal A is then given by 

pfor_n(l) f^a) _ 19^(1) «^ 

— °^AAAV ^^AAAA — ^'^^AV 

The interaction e^^y is deduced from the experimental value of the vacancy forma- 
tion energy in pure Al. For the interaction e^rv assume that the vacancy formation 
energy in the fee structure is the same as in the hep one, these two structures being 
quite similar. The only experimental information we have concerning this energy is 
Ey^ > 1.5 eV [30]. We thus use the ab-initio value calculated by Le Bacq et al. [31], 
Ey^ = 2.07 eV. This value is calculated at the equilibrium volume of Zr and cannot be 
used directly to obtain e^^y as this interaction should correspond in our model to the 
equilibrium lattice parameter of pure Al. We have to use instead the vacancy formation 
enthalpy 

hI°' = El;' + P5VLi;\ (14) 

where 6Vt{;'' = -1.164 is the vacancy formation volume in pure Zr [30], and P is 
the pressure to impose to Zr to obtain a lattice parameter equal to the one of Al. P is 
calculated from the bulk modulus i? = 91 GPa of fee Zr and the equilibrium volumes 
of Al and Zr, = 16.53 and = 23.36 A^, these three quantities being 
obtained from our FP-LMTO calculations. This gives us the value = 1.88 eV for 
the vacancy formation enthalpy in pure Zr at the lattice parameter of pure Al. 

We use the experimental value of the divacancy binding energy i?2y = 0-2 eV [30] 
in order to compute a vacancy- vacancy interaction, eyy = i^gy — ^aiai + "^^Aiv 
do not include this interaction and set it equal to zero instead, we obtain the wrong 
sign for the divacancy binding energy, divacancies being thus more stable than two 
monovacancies. This does not affect our Monte Carlo simulations as we only include 
one vacancy in the simulation box, but this will have an influence if we want to build 
a mean field approximation of our diffusion model. 

We thus managed to add vacancy contributions to our thermodynamic description 
of the Al-Zr binary system. Using the breakdown [3 of the first nearest neighbor tetrahe- 
dron interaction, we can obtain the effective energies corresponding to the 15 different 
configurations a tetrahedron can have in the ternary system. These effective energies 
are presented in table lU for the cases where energy correction due to order on first 
nearest neighbor triangle and tetrahedron are assumed different from zero or not (J3 
and J4 given by table 01 or J3 = J4 = 0). 



(13) 
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Table 4: Effective energies of tlie first nearest neiglibor tetraliedron for Al-Zr-V ternary 
system. Tlie set witli order correction corresponds to tlie values J3 and J4 given by the 
cluster expansion of table El and the set without order correction assumes J3 = J4 = 0. 



Configuration 


Effective 


energy (eV) 




with order 


without order 




correction 


correction 


Al Al Al Al 


-1.680 


-1.680 


Al Al Al Zr 


-2.257 


-2.214 


Al Al Zr Zr 


-2.554 


-2.554 


Al Zr Zr Zr 


-2.707 


-2.698 


Zr Zr Zr Zr 


-2.647 


-2.647 


Al Al Al V 


-1.174 


-1.174 


Al Al Zr V 


-1.567 


-1.561 


Al Zr Zr V 


-1.748 


-1.754 


Zr Zr Zr V 


-1.751 


-1.751 


Al Al V V 


-0.518 


-0.518 


Al Zr V V 


-0.758 


-0.758 


Zr Zr V V 


-0.804 


-0.804 


Al V V V 


+0.288 


+0.288 


Zr V V V 


+0.194 


+0.194 


WW 


+ 1.243 


+1.243 



With this set of thermodynamic parameters, we calculate the binding energy be- 
tween a Zr solute atom and a vacancy in Al, 

pbm _ 9 / (1) , Al) _ ii) _ fi) A 

^ZrV — ^ [^AIAIAIAI ^AlAlZrV ^AlAlAlZr ^AIAIAIV J ■ \^'~') 

The value obtained, = +0.369 eV, agrees with the experimental observation that 

there is no attraction between Zr solute atoms and vacancies in Al [30,32]. 



3.4. Migration barriers 

Diffusion occurs via vacancy jumps towards one of its twelve first nearest neighbors. 
The vacancy exchange frequency with a neighbor of type A {A =A1 or Zr) is given by 

~~k^)' ^^^^ 

where va is an attempt frequency and the activation energy is the energy change 
required to move the A atom from its initial stable position to the saddle point position. 
It is computed as the difference between the contribution e^* of the jumping atom to the 
saddle point energy and the contributions of the vacancy and of the jumping atom to the 
initial energy of the stable position. This last contribution is obtained by considering 
all bonds which are broken by the jump, i.e. all pair interactions the vacancy and the 
jumping atoms are forming as well as all order corrections on triangles and tetrahedrons 
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containing the jumping atom, 

j j^V jk jkl j 

The attempt frequency pa and the contribution e^of the jumping atom to the saddle 
point energy can depend on the configuration [2] . Nevertheless, we do not have enough 
information to see if such a dependence holds in the case of Al-Zr alloys. We thus 
assume that these parameters depend only on the nature of the jumping atom, which 
gives us four purely kinetic parameters to adjust. 

The contribution of Al to the saddle point energy, e^^, is deduced from the exper- 
imental value of the vacancy migration energy in pure Al, Ey^^ = 0.61 eV [30], and 
the attempt frequency u^i from the experimental Al self-diffusion coefficient, D^i* = 
DoBxp {—Q/ksT), the self-diffusion activation energy Q being the sum of the vacancy 
formation and migration energies in pure Al and Dq = 1.73 x 10~^ m^.s~^ [33]. 

To calculate uzr and e^^, we use the experimental value^ of the diffusion coefficient of 
Zr impurity in Al, Dzr* = 728 x lO'^exp (-2.51 eV/kBT) m^.s^^ [33,34]. The kinetic 
parameters can be deduced from this experimental data by using the five frequency 
model for solute diffusion in fee lattices [35], if we make the assumption that there is no 
correlation effect. We check afterwards that such an assumption is valid: at T = 500 K 
the correlation factor is fzr* = 1 and at T = 1000 K fzr* = 0.875. Correlation effects 
are thus becoming more important at higher temperature but they can be neglected in 
the range of temperature used in the fitting procedure. 

Table 5: Kinetic parameters for a thermodynamic description of Al-Zr binary with 
and without energy corrections due to order on first nearest neighbor triangle and 
tetrahedron. 





with order 




without order 






correction 




correction 






-8.219 eV 




-8.219 eV 




p^P 
^Zr 


-11.286 eV 




-10.942 eV 






1.36 X 10^^ 


Hz 


1.36 X 10^^ 


Hz 


yzr 


4.48 X 10^^ 


Hz 


4.48 X 10^^ 


Hz 



So as to study the infiuence on kinetics of energy corrections due to order on triangles 
and tetrahedrons, we fit another set of kinetic parameters corresponding to a thermo- 
dynamic description of Al-Zr binary with only pair interactions [i.e. J3 = J4 = 0, or 
equivalently e^j], = e^j^^ = 0). This other set of kinetic parameters presented in tabled] 
reproduces as well coefficients for Al self-diffusion and for Zr impurity diffusion, the only 
difference being that these kinetic parameters correspond to a simpler thermodynamic 
description of Al-Zr binary. 



Diffusion coefficient measured in the temperature range between 800 and 910 K. 
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4. DIFFUSION IN SOLID SOLUTION 

Diffusion in Al-Zr solid solutions can be fully characterized by the tracer correlation 
coefficients /ai and fzr and by the phenomenological Onsager coefficients Laiai, Laizt, 
and LzrZr- These coefficients link ffuxes of diffusing species, Jai and Jzr, to their 
chemical potential gradients [36, 37] through the relations 

Jai = -LaiaNij^'ai/^bT - LAlZryiJ^Zr/^BT 

Jzr = -LaIZt^I^'ai/^bT - LzrZr^f^'zr/^BT. (18) 

Chemical potentials entering these equations are relative to the vacancy chemical po- 
tential, fi'j^i = II Ai — l^v and ii'zr = ^^Zr — fi'V- We use to express diffusion fluxes the 
Onsager reciprocity condition, Laizv = LzrAi- 

These coefficients can be used in finite-difference diffusion code so as to study "in- 
dustrial" processes where diffusion is involved (precipitation, solidification, homogeniza- 
tion, . . . ) [38]. One way to obtain these coefficients is to adapt Calphad methodology 
to kinetics, i.e. to guess an expression for Lab describing its variation with temper- 
ature and composition of the alloy and to adjust the model parameters on a large 
kinetic database [39,40]. On the other hand one can use an atomistic model as the 
one described in chap. ID to obtain the phenomenological coefficients [8,9,41]. If one 
carefully applies the same mean field approximation for thermodynamics and kinetics 
it is possible to get the whole Onsager matrix and not only diagonal terms and to catch 
all correlation effects [42]. Such an approach compared to the previous one does not 
need a huge experimental database. Moreover, as it is based on a realistic description 
of diffusion at the atomic scale, it appears safer to extrapolate kinetic quantities out of 
the range (composition or temperature) used in the fitting procedure. 

In this study, we do not use any mean -field approximation to calculate phenomeno- 
logical coefficients, but obtain them directly from kinetic Monte Carlo simulations by 
using generalization of the Einstein formula for tracer diffusion due to Allnatt [37,43] 

L,,= <^5£_^, A.B = Al,Zr. (19) 

where the brackets indicate a thermodynamic ensemble average and AR^i is the sum 
of total displacement Ar^ of all atoms i of type A during time At, 

AR^ = ^Ar,. (20) 

We use residence time algorithm to run kinetic Monte Carlo calculations. The 
simulation box contains 125000 lattice sites, one of this site being occupied by a vacancy. 
Sum of total displacements AR^^ and AR^^ in equation UHl are computed for a time 
interval corresponding to ~ 10^ vacancy jumps, and their thermodynamic averages are 
obtained through simulations of 10^ vacancy jumps. Such a big number of jumps is 
necessary to converge thermodynamic averages entering in the calculation of Laizt and 
LzrZr, whereas Laiai converges more quickly. This is due to the difference of diffusion 
coefficients between Al and Zr. 

Results of calculations are presented in figure |21 for two different temperatures, 
T = 1000 K and T = 900 K, and different Zr concentration from to 8 at.%. For the 
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Figure 2: Onsager coefficients Laiai and LzrZr- Squares and solid lines correspond to T = 1000 K 
and circles and dashed lines to T = 900 K. The vertical lines indicate the corresponding solubility 
limit obtained from CVM calculations. Full symbols correspond to the set of parameters with order 
corrections on triangles and tetrahedrons and open symbols to the set without order corrections. 

off-diagonal coefficient Laizt of Onsager matrix, dispersion is too important to get a 
precise value of thermodynamic average^. We interpret this as an indication that this 
coefficient can be neglected in this range of temperature and concentration. 

Onsager coefficients are calculated for Zr concentration corresponding to the stable 
as well as to the metastable solid solution, the limit being given by the CVM calculations 
of chap.|2I2Il For calculations in the metastable solid solution, thermodynamic averages 
are computed during the incubation stage of precipitation kinetics when no stable 
precipitate is present in the simulation box (chap. EU- Laiai behavior deviates only 
slightly from its linear extrapolation from the stable solid solution, but for LzrZr it 



"^Laizt = ± 10^12 ni2.s-i at T = 1000 K and Laizt = ± 10"" m^.s'^ at T = 900 K 
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seems that no extrapolation from the stable to the metastable solid solution is possible. 

So as to see the influence of triangles and tetrahedrons interactions, we ran simu- 
lations with only pair interactions considering the corresponding kinetic parameters of 
table El One can directly see on figure |21 that these two sets of parameters reproduce 
the same experimental data, i.e. the self-diffusion coefficient 

Dai' = /o lim Laiai, (21) 
where /o = 0.78145 for a fee lattice, and the Zr impurity diffusion coefficient 

Dzr* = hm LzrZr/Czr- (22) 

Czr^O 

Order corrections mainly affect LzrZr- This coefficient is slightly lower when one 
considers energy corrections due to order on triangles and tetrahedrons. The differ- 
ence increases with Zr concentration and thus in the metastable solid solution: these 
order corrections lead to a slight slowdown of Zr diffusion. The two thermodynamic 
models are equivalent at these temperatures (c/. phase diagram on Fig. [U (a)). As a 
consequence kinetic behaviors obtained from them are really close. 



1 .6e-07 




2 4 6 8 10 12 

%Zr 

Figure 3: Onsager coefficients Laiai, LzrZr, and Laizt calculated at T = 3000 K. Full symbols and 
solid lines correspond to the set of parameters with order corrections on triangles and tetrahedrons 
and open symbols and dashed lines to the set without order corrections. 

At higher temperatures, triangle and tetrahedron interactions change the phase di- 
agram (Fig. n (a)). This thermodynamic influence leads to a kinetic change too: at 
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T = 3000 K, Onsager coefficients are lower when considering these multisite interac- 
tions (Fig. 121). One should notice that correlation effects cannot be neglected at this 
temperature as L^izr is far from being null. Thus one is not allowed anymore to as- 
sume Onsager matrix as diagonal. With triangle and tetrahedron interactions, AlsZr 
precipitate is more stable, which means that order effects are stronger. The kinetic 
corrolary of this thermodynamic influence is that they slow down diffusion. 

5. KINETICS OF PRECIPITATION 

Precipitation kinetics have been obtained by Monte Carlo simulations for four differ- 
ent supersaturations of the solid solution {C^zr — 5; 6, 7, and 8 at.%) at T = 1000 K. At 
this temperature, the equilibrium concentration is C^^r ~ 2-1 at.%. The simulation box 
contains 125000 lattice sites and its starting configuration is a completely disordered 
(random) solid solution. 

5.1. Short range order parameters 



0.3 




0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 0.6 

t (ms) t (ms) 

Figure 4: Evolution of second nearest neighbor short range order of Zr atoms, a\^, at T = 1000 K 
and four different nominal concentrations C^^. Full and dotted lines are respectively for the set of 
parameters with and without order corrections on triangles and tetrahedrons. 

The quantities of interest to follow the global evolution of precipitation during the 
simulation are Warren- Cowley short range order (SRO) parameters [5]. SRO parame- 
ters for first-nearest neighbors evolve too quickly to give any really significant informa- 
tion on precipitation state. During simulation ffist steps, Zr atoms surround themselves 
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with Al. Once this local equilibrium for first nearest neighborhood is reached, the corre- 
sponding SRO parameters do not evolve anymore. On the other hand, SRO parameters 
for second nearest neighbors slowly evolve until the end of the simulation. For Zr atoms, 
it is defined as 

2 {'Pn'')zr,2 ~ ^Zr 

^zr = ^3^1 , (23) 

where {Vn^) stands for the average of occupation numbers pf'' on all second nearest 
neighbors of Zr atoms. For a randomly distributed configuration of the alloy (initial 
configuration) = 0, whereas for the LI2 structure = 1. Looking at fig. IH one 
sees that evolves more quickly with the set of parameters with only pair interactions 
than with triangle and tetrahedron interactions. At first glance, this is in agreement 
with the slight difference on LzrZr measured in the metastable solid solution at this 
temperature (Fig. |2)) for the two set of parameters. So as to see if the difference of 
precipitation kinetics can be understood only in terms of a difference of diffusion speed 
or is due to another factor, we measure the nucleation rate in our simulations and 
interpret it with classical theory of nucleation [1,44,45]. 

5.2. Precipitate critical size 

We first need to give us a criterion to decide which atoms are belonging to LI2 
precipitates. As stable precipitates are almost perfectly stoichiometric at T = 1000 K 
(chap. l2.H.j) . we only look at Zr atoms and consider for each Zr atom in LI2 precipitate 
that three Al atoms are belonging to the same precipitate. Zr atoms are counted as 
belonging to LI2 precipitates if all their twelve first nearest neighbors are Al atoms and 
at least half of their six nearest neighbors are Zr atoms. Moreover, we impose that at 
least one Zr atom in a precipitate has its six second nearest neighbors being Zr, i. e. has 
a first and second nearest neighborhood in perfect agreement with the LI2 structure. 

Classical theory of nucleation predicts there is a critical radius, or equivalently a 
critical number i* of atoms, below which precipitates are unstable and will re-disolve 
into the solid solution and above which precipitates will grow, i* is obtained by con- 
sidering the competition between the interface free energy a and the nucleation free 
energy per atom AG", 

2tx ( a^a 



AG" 



(24) 



Clusters of size i < i* are considered to be local variations of the solid solution compo- 
sition and thus are not counted as LI2 precipitates. 

5.3. Nucleation free energy 

The nucleation free energy per atom entering equation |^ is given by [44,45] 

AG" = I i^^AiiCZ) - miCl)) + \ if^zriCZ) - /iz.(Gi)) , (25) 

where fiAi{Czr) and /xzr(Czr) are the chemical potentials of respectively Al and Zr com- 
ponents in the solid solution of concentration Czr, C^zr equilibrium concentration 
of the solid solution, and G^^ the nominal concentration. The factors 3/4 and 1/4 
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Figure 5: Nucleation free energy AG" at T = 1000 K for different concentration of the solid solution. 
Square symbols correspond to CVM-TO calculation and the line to the ideal solution approximation 
(eq. I26|l . Full and open symbols are respectively for the set of parameters with and without order 
corrections on triangles and tetrahedrons. 

arises from the stoichiometry of the precipitating phase AlaZr. Usually the nucleation 
free energy is approximated by 

AG- = 4^BTlog^-^ + ^^BTlog^ (26) 

which is obtained by considering in equation 1221 the expressions of the chemical poten- 
tials for an ideal solution. As at T = 1000 K we obtained the same solubility limit, 
= 2.1 at.%, with or without triangle and tetrahedron interactions (c/. phase dia- 
gram on fig. Q), the approximation 1^ cannot be used to see if these interactions have 
any influence on the nucleation free energy. Therefore we use CVM-TO to calculate 
chemical potentials entering expression 123 Looking at figure El one should notice that 
the ideal solution approximation would have lead to an overestimation of AG", the 
error being ~ 10% for the maximal supersaturation considered. With CVM-TO, we 
do not obtain any change in the value of the nucleation free energy depending we are 
considering or not order corrections for first nearest neighbor triangle and tetrahedron. 
Thus slowdown of precipitation kinetics with these corrections cannot be explained by 
a decreasing of the nucleation free energy. 

5.4. Interface free energy 

To determine the precipitate critical size i* using expression we need to know 
the value of the interface free energy a too. We calculate this energy at K for 
different orientations of the interface. We therefore do not consider any configurational 
entropy and simply obtain the interface energy by counting the number by area unit 
of wrong "bonds" compared to pure Al and AlsZr in LI2 structure. For (100) and 
(110) interfaces there is an ambiguity in calculating such an energy as two different 
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planes, one pure Al and the other one of stoichiometry Ali/2Zri/2, can be considered as 
interface. Considering LI2 precipitates as stoichiometric will guarantee that to any type 
of the two possible interfaces is associated a parallel interface of the other type. Thus 
for (100) and (110) interfaces, we consider the average of these two different interface 
energies to be meaningful for the parameter a entering in classical theory of nucleation. 
For (111) interface, as only one interface of stoichiometry Al3/4Zri/4 is possible, we do 
not obtain such an ambiguity. The energies corresponding to these different interfaces 
are 

1 1 iS^') _ _ ^(2) 

cTioo = ^^110 = ^111 = 2a2^ ~' "^"^^^^ " ^^'^ 

These interface energies only depend on second nearest neighbor interactions and there- 
fore are the same with or without order corrections on first nearest neighbor triangle 
and tetrahedron. To determine the critical size of precipitates with equation |^ we use 
an interface free energy slightly higher than aioo, o^cr = 64.1 meV. With this interface 
free energy, nucleation rate obtained from Monte Carlo simulations are in better agree- 
ment than with aioo (Fig- IZ|)- As precipitates observed in Monte Carlo simulations do 
not exhibit sharp interfaces, this is quite natural to have to use an energy higher than 
the minimal calculated one. 



5.5. Nucleation rate 

Critical size for precipitates obtained from these nucleation and interface free ener- 
gies are respectively i* = 187, 104, 76, and 57 atoms for the different nominal concen- 
trations C^r = 5, 6, 7, and 8 at.%. We use these critical sizes to determine the number 
Np of supercritical precipitates contained in the simulation boxes, their average size 
as well as the concentration of the solid solution Czr- The variation with time of 
these quantities are shown on fig. IHl for the simulation box of nominal concentration 
= 8 at.%. After an incubation time, one observes a nucleation stage where the 
number of precipitates increases linearly until it reaches a maximum. We then enter 
into the growth stage: the number of precipitates does not vary and their size is in- 
creasing. At last, during the coarsening stage, precipitates are still growing but their 
number is decreasing. For this concentration, one clearly sees that precipitation kinet- 
ics is faster with only pair interactions as the number of precipitates is increasing more 
rapidly. Moreover precipitates have a bigger size than with triangle and tetrahedron 
order corrections. 

The steady-state nucleation rate J'^* is measured during the nucleation stage, when 
the number of precipitates is varying quite linearly with time. Slowdown of precipitation 
kinetics with triangle and tetrahedron order corrections can be seen on the steady-state 
nucleation rate (Fig. [7j): without these corrections J'^* is about two times higher than 
when these corrections are included. 

In classical theory of nucleation, the steady-state nucleation rate is given by the 
expression [44], 

AC* 

r^ = NoZp*exp-^, (27) 

where A'^o is the number of nucleation sites, i.e. the number of lattice sites (A'^o = 125000 
for Monte Carlo simulations), AG* is the nucleation barrier and corresponds to the free 
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Figure 6: Kinetics of precipitation for a nominal Zr concentration — 8 at.%: evolution with 
time of the number iVp of precipitates in the simulation box, of precipitates average size < i >p, and of 
Zr concentration in the solid solution. Full and dotted lines are respectively for the set of parameters 
with and without order corrections on triangles and tetrahedrons. 



energy of a precipitate of critical size i* 



AG* = V, (2^ 

Z is the Zeldovitch factor and describes size fluctuations of precipitates around i*, 



1 AG"2 



27r {a^ay/^VkT 



(29) 



and (3* is the condensation rate for clusters of critical size i*. Assuming the limiting 
step of the adsorption is the long range diffusion of Zr atoms in the solid solution, the 
condensation rate is [44] 

Zr diffusion coefficient is obtained from our measure of Onsager coefficients in the 
metastable solid solution (chap. Assuming that vacancies are at equilibrium (/iy = 
0), its expression is [35-37] 

n -(r r \ 1 ^^^r rsn 



CI J ksTdCl 
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Figure 7: Evolution of the steady-state nucleation rate J** with the nominal concentration for 
T = 1000 K. Full and open symbols are respectively for the set of parameters with and without order 
corrections on triangles and tetrahedrons. The full line corresponds to the nucleation rate predicted by 
classical theory of nucleation with a?(j — 64.1 meV and the dotted line with o?a = a^cioo = 57.0 meV. 
J** is normalized by the number of lattice sites in the simulation box, = 125000. 



We obtain the thermodynamic factor dfizr/dC^r using CVM-TO calculations. This 
factor is the same with or without order corrections on triangles and tetrahedrons. 
Therefore, the only difference these corrections induce on Zr diffusion arises from LzrZr- 
In classical theory of nucleation, the diffusion coefficient entering in the expression EDI 
of the condensation rate is only a scaling factor for time and does not have any other 
influence on kinetics. As a consequence the steady-state nucleation rate varies linearly 
with Zr diffusion coefficient as it clearly appears when combining equations IHUl and 1271 
Thus small variations of LzrZr with the set of parameters used do not allow to explain 
the difference of the nucleation rate: with order corrections, LzrZr is far from being 
half the value it is with only pair interactions (Fig. Ej). Thus slowdown of precipitation 
kinetics is not due to a slowdown of Zr diffusion. 

One possible explanation would be a difference of the interface free energy a. J** 
is really sensitive to this parameter and one only needs a small decrease of a to obtain 
a higher nucleation rate (see on fig. [7|the decrease of J'^* when a^a is going from 57.0 
to 64.1 meV). Such a decrease would explain too why precipitates have a bigger size 
with only pair interactions. At T = K, we obtain the same interface energy for all 
directions considered with the two sets of parameters, but at finite temperature the 
configurational entropy could lead to a difference of interface free energy. Nevertheless, 
this needs to be confirmed, by CVM calculations or using the Cahn-Hilliard method [46] 
for instance. Another possible explanation to understand the different kinetic pathways 
would be a different mobility of small clusters with the two sets of parameters. But 
this would be quite surprising, as we do not expect small clusters to be really mobile 
because of the repulsion between vacancy and Zr solute atoms. 
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6. CONCLUSIONS 

We built an atomistic kinetic model for Al-Zr binary system using ab-initio calcu- 
lations as well as experimental data. So as to be as realistic as it should be at this 
atomic scale, this model describes diffusion through vacancy jumps. Thanks to ab- 
initio calculations we could improve usual thermodynamic descriptions based on pair 
interactions and incorporate multisite interactions for clusters containing more than 
two lattice points so as to consider dependence of bonds with their local environment. 

At temperatures lower than 1000 K, these energetic corrections due to local order 
do not modify thermodynamics: the phase diagram does not change when one does not 
consider these order corrections. For higher temperatures they lead to a stabilization 
of the ordered structure LI2. 

Concerning diffusion in the solid solution, these order corrections on first nearest 
neighbor triangle and tetrahedron do not really change the Onsager matrix, and thus 
diffusion characteristics. They just lead to a slight slowdown of Zr diffusion in the 
metastable solid solution. When looking at higher temperatures, the slowdown of Zr 
diffusion is more important. 

For precipitation, kinetics are slower with these interactions. The slowdown is too 
important to be related to the small decrease of Zr diffusion in the metastable sohd 
solution at the same temperature. One possibility would be a change of configurational 
entropy contribution to interface free energy. 
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A DETAILS OF AB-INITIO CALCULATIONS 

Ab initio calculations were carried out using a full-potential linear-muffin-tin-orbital 
(FP-LMTO) method [17-19] in the version developed by Methfessel and Van Schilf- 
gaarde [47]. The basis used contained 22 energy independent muffin-tin-orbitals (MTO) 
per Al and Zr site: three k. values for the orbitals s and p and two k, values for the or- 
bitals d where the corresponding kinetic energies were = 0.01 Ry (spd), 1.0 Ry (spd), 
and 2.3 Ry (sp). A second panel with a basis composed of 22 energy independent MTO 
with the same kinetic energies was used to make a correct treatment of the 4p semicore 
states of Zr. The same uniform mesh of points was used to make the integrations in 
the Brillouin zone for valence and semicore states. The radii of the muffin-tin spheres 
were chosen to have a compactness of 47.6% for Al sites and 54.1% for Zr sites. Inside 
the muffin-tin spheres, the potential is expanded in spherical harmonics up to / = 6 
and in the interstitial region spherical Hankel functions of kinetic energies = 1 Ry 
and 3.0 Ry were fitted up to i = 6. The calculations were performed in the generalized 
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gradient approximation (GGA) [48,49] and the parameterization used was the one of 
Perdew et al. [50]. 
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